Modèle d’ensembles échangeable
Un Exemple de Températures pour l’approche normale
Eric Parent
17 juillet 2019
1 Introduction
Jacques et moi avons beaucoup communiqué sur les modèles d’ensemble échangeables, mais absolument rien valorisé en terme d’article vis à vis de la communauté d’hydrologie statistique que nous nous plaisons tant à critiquer. D’autre part, nos accords (ainsi que désaccords) théoriques sur l’intérêt du caractère échangeable pour modéliser les ensembles et du Bayesian Forecasting System de Krzysztofowicz afin de construire une prédictive resteront du domaine des plaisantes joutes epistolaires limitées à deux crocodiles bayésiens si nous ne montrons pas leur caractère opératoire sur de “vraies” jeux de données. Je propose donc de traiter de tels jeux de données pour lesquels on mettra en évidence les améliorations concrètes et pratiques que ces méthodes apportent (au delà de leur caractère théorique rationel et cohérent). C’est l’objectif de ce papier, écrit en Markdown html notebook dans l’intention de donner une illustration informatique et applicative en même temps que du texte plus théorique. J’imagine en fait deux papiers de ce type:
le premier ici sur le modèle normal échangeable inspiré du modèle à deux pivots (voir pour le jus théorique le texte de Jacques ENSEMBLESAISON.tex) que l’on peut illustrer sur des données de températures d’EDF-DTG,
le second, à venir, sur des données zéros-inflatées de précipitations, où l’on mettrait en oeuvre les modèles Bernoulli-Gamma et lois des fuites. Ces papiers pourront servir de base à des articles publiables à terme… que j’espère pas très lointain.
2 Analyses Descriptives
2.1 Les données disponibles
rm(list=ls())
# setwd("~/Documents/courbariaux/WORKLUC/TempEnsemblesNormal")
load("Ain@Vouglans.Rdata")
library(tidyverse)Registered S3 method overwritten by 'dplyr':
method from
print.rowwise_df
Registered S3 methods overwritten by 'ggplot2':
method from
[.quosures rlang
c.quosures rlang
print.quosures rlang
[37m── [1mAttaching packages[22m ─────────────────────────────────────── tidyverse 1.2.1 ──[39m
[37m[32m✔[37m [34mggplot2[37m 3.1.1 [32m✔[37m [34mpurrr [37m 0.3.2
[32m✔[37m [34mtibble [37m 2.1.3 [32m✔[37m [34mdplyr [37m 0.8.1
[32m✔[37m [34mtidyr [37m 0.8.3 [32m✔[37m [34mstringr[37m 1.4.0
[32m✔[37m [34mreadr [37m 1.3.1 [32m✔[37m [34mforcats[37m 0.4.0[39m
[37m── [1mConflicts[22m ────────────────────────────────────────── tidyverse_conflicts() ──
[31m✖[37m [34mdplyr[37m::[32mfilter()[37m masks [34mstats[37m::filter()
[31m✖[37m [34mdplyr[37m::[32mlag()[37m masks [34mstats[37m::lag()[39mlibrary(tidyselect)
library(lubridate)
Attaching package: ‘lubridate’
The following object is masked from ‘package:base’:
datelibrary(R2jags)Loading required package: rjags
Loading required package: coda
Linked to JAGS 4.3.0
Loaded modules: basemod,bugs
Attaching package: ‘R2jags’
The following object is masked from ‘package:coda’:
traceplotlibrary(verification)Loading required package: fields
Loading required package: spam
Loading required package: dotCall64
Loading required package: grid
Spam version 2.2-2 (2019-03-07) is loaded.
Type 'help( Spam)' or 'demo( spam)' for a short introduction
and overview of this package.
Help for individual functions is also obtained by adding the
suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
Attaching package: ‘spam’
The following objects are masked from ‘package:base’:
backsolve, forwardsolve
Loading required package: maps
Attaching package: ‘maps’
The following object is masked from ‘package:purrr’:
map
See https://github.com/NCAR/Fields for
an extensive vignette, other supplements and source code
Loading required package: boot
Loading required package: CircStats
Loading required package: MASS
Attaching package: ‘MASS’
The following object is masked from ‘package:dplyr’:
select
Loading required package: dtw
Loading required package: proxy
Attaching package: ‘proxy’
The following object is masked from ‘package:spam’:
as.matrix
The following objects are masked from ‘package:stats’:
as.dist, dist
The following object is masked from ‘package:base’:
as.matrix
Loaded dtw v1.20-1. See ?dtw for help, citation("dtw") for use in publication.library(GGally)Registered S3 method overwritten by 'GGally':
method from
+.gg ggplot2
Attaching package: ‘GGally’
The following object is masked from ‘package:dplyr’:
nasaSuite à la thèse de Marie Courbariaux, je dispose des données journalières de températures fournies par EDF-DTG aux stations suivantes:
l’ Ain à Vouglans,
le Buech à Chambons,
le Drac à Sautet.
Je propose de travailler dans toute la suite sur l’Ain à Vouglans, mais les programmes développés sont directement transposables aux deux autres stations. Nous disposons de l’historique journalier sur la période 1953 à 2015 dans le data.frame dhisto, ainsi que des \(50\) membres du Centre Européen de Prévision pour les années 2005 à 2008 (4ans) puis de 2O11 à 1015 (5 ans) pour des prévisions jusqu’à une échéance de 9 jours dans le data.frame dtout.
dhisto %>% glimpseObservations: 23,010
Variables: 7
$ Date [3m[38;5;246m<date>[39m[23m 1953-01-01, 1953-01-02, 1953-01-03, 1953-01-04, 1953-01-0…
$ Qobs [3m[38;5;246m<dbl>[39m[23m 33.0, 30.0, 29.1, 24.7, 23.7, 23.9, 22.5, 21.9, 19.6, 18.3…
$ PS [3m[38;5;246m<dbl>[39m[23m 0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 1.1, 0.0, 4.7, 3.9, 0.0, 0.0…
$ T [3m[38;5;246m<dbl>[39m[23m -0.8, 0.0, -1.4, -1.9, -1.4, -1.8, -3.9, -2.6, -1.6, -1.6,…
$ Year [3m[38;5;246m<dbl>[39m[23m 1953, 1953, 1953, 1953, 1953, 1953, 1953, 1953, 1953, 1953…
$ Month [3m[38;5;246m<dbl>[39m[23m 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ Calendaire [3m[38;5;246m<dbl>[39m[23m 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,…range(dhisto$Date)[1] "1953-01-01" "2015-12-31"dtout %>% head()unique(dtout$Year) [1] 2005 2006 2007 2008 2009 2011 2012 2013 2014 2015 2016unique(dtout$Echeance)[1] 1 2 3 4 5 6 7 8 92.2 Analyse climatologique
On peut illustrer les fluctuations saisonnières des températures journalières au cours de l’année:
dhisto %>% mutate(an= as_factor(Year)) %>% ggplot(aes(x = Calendaire, y= T, color = an)) + geom_line(alpha=0.5)
horsain = ymd("2003-10-15")Ceci met en évidence un comportement très sûrement aberrant en date du 2003-10-15, horsain que l’on enlève
dhisto %>% filter(Date!= horsain) %>% mutate(an= as_factor(Year)) %>% ggplot(aes(x = Calendaire, y= T)) + geom_point(aes( color = an),alpha=0.5, shape ='.', show.legend = FALSE)+geom_smooth()+labs(x='Jours Calendaires', y='Températures', title="Climatologie à Vouglans")
On va alors construire une température moyenne de réference et son écart-type associé, éventuellement lissés pour chaque jour calendaire de l’année. Pour construire ces estimateurs périodiques, on régresse sure une base des 4 premières harmoniques la température interannuelle moyenne et le logarithme de l’écart-type journalier de la température.
dhisto %>% filter(Date!= horsain, Calendaire != 366) %>% mutate(an= as_factor(Year)) %>% group_by(Calendaire) %>% summarize (T_moy_ref= mean(T),T_std_ref= var(T)^0.5 ) -> T_ref
modlin<-lm(T_moy_ref~ 1+ sin(2*pi*Calendaire/365)+ cos(2*pi*Calendaire/365)+
sin(4*pi*Calendaire/365) +
cos(4*pi*Calendaire/365)+
sin(6*pi*Calendaire/365) +
cos(6*pi*Calendaire/365)+
sin(8*pi*Calendaire/365)+
cos(8*pi*Calendaire/365), data=T_ref
)
modlinStd <-lm(log(T_std_ref)~ 1+ sin(2*pi*Calendaire/365)+ cos(2*pi*Calendaire/365)+
sin(4*pi*Calendaire/365) +
cos(4*pi*Calendaire/365)+
sin(6*pi*Calendaire/365) +
cos(6*pi*Calendaire/365)+
sin(8*pi*Calendaire/365)+
cos(8*pi*Calendaire/365), data=T_ref
)
T_ref %>% mutate(T_moy_ref_lis=modlin$fitted.values,
T_std_ref_lis=exp(modlinStd$fitted.values) )-> T_ref
T_ref%>%
ggplot(aes(x=Calendaire))+geom_line(aes(y=T_moy_ref_lis)) + geom_point(aes(y=T_moy_ref))+labs(caption="Moyenne calendaire (référence climatologique)")
T_ref %>%
ggplot(aes(x=Calendaire))+geom_line(aes(y=T_std_ref_lis)) + geom_point(aes(y=T_std_ref))+labs(caption="Ecart-Type calendaire (référence climatologique) ")
2.3 Analyse descriptive de l’ensemble du centre européen de prévision
L’analyse climatologique précédente permet de centrer et réduire les températures pour travailler en anomalies journalières. Ces anomalies sans dimension sont alors considérées comme gaussiennes, marginalement \(N(0,1)\). On va laisser tomber les membres individuels pour en conserver les statistiques exhaustives, moyenne et variance empiriques dans le dataframe de travail d. Construisons d’abord une base de données (en ajoutant la climatologie).
dtout %>% mutate(Xbar=rowMeans(dplyr::select(.,starts_with("Run"))),
V2=apply(dplyr::select(.,starts_with("Run")),1,var)) %>%
dplyr::select(-starts_with("Run")) ->d2.3.1 Relation entre moyenne de l’ensemble et température observée
Représentons en premier lieu le lien entre température observée et température moyenne prévue à diverses échéances (après centrage par la moyenne climatologique calendaire et réduction par l’écart-type climatologique calendaire )
d %>% left_join(T_ref %>% dplyr::select(-T_std_ref,-T_moy_ref)) ->dJoining, by = "Calendaire" d %>% ggplot(aes(x=(Obs-T_moy_ref_lis)/T_std_ref_lis,y=(Xbar-T_moy_ref_lis)/T_std_ref_lis,color=as.factor(Month)) )+geom_point()+geom_abline(slope=1,intercept=0) +facet_wrap(~Echeance)+labs(x="Anomalie de température", y="Anomalie de prévision", color="Mois")
Pour la moyenne d’ensemble, on constate que le biais et la dispersion vis à vis de la cible augmentent avec l’échéance, comme on s’y attend.
d %>% mutate(x=(Obs-T_moy_ref_lis)/T_std_ref_lis,y=(Xbar-T_moy_ref_lis)/T_std_ref_lis) %>% group_by(Echeance) %>% summarise(biais=mean(x-y, na.rm=T), variance=var(x-y, na.rm=T), sce=mean((x-y)^2, na.rm=T))Pour une échéance donnée, il y a également des fluctuations selon le mois considéré, mais ces différences restent relativement raisonnables. A titre d’exemple, faisons le calcul pour une échéance de 4 jours:
d %>% filter(Echeance ==4) %>% mutate(x=(Obs-T_moy_ref_lis)/T_std_ref_lis,y=(Xbar-T_moy_ref_lis)/T_std_ref_lis) %>%
group_by(Month) %>% summarise(biais=mean(x-y, na.rm=T), variance=var(x-y, na.rm=T),sce=mean((x-y)^2, na.rm=T))2.3.2 Relation entre variance de l’ensemble et température observée
Par contre, la variance inter-membre ne semble guère reliée à la température observée
d %>% ggplot(aes(x=(Obs-T_moy_ref_lis)/T_std_ref_lis,y=log(V2/(T_std_ref_lis^2)),color=as.factor(Month)) )+geom_point()+facet_wrap(~Echeance)+labs(x="anomalie de température", y="Variance relative (enLog)", color="Mois")
Comme le montre la figure précédente, la valeur moyenne (géométrique) de la variance inter-membre augmente avec l’échéance jusqu’à atteindre 42% de la variance climatologique à l’échéance 9 jours, comme on s’y attend, mais pas la dispersion de son logarithme, qui, au contraire, diminue .
d %>% mutate(y=log(V2/(T_std_ref_lis^2))) %>% group_by(Echeance) %>% summarise(moygeomvariance=exp(mean(y, na.rm=T)), VarLogVariance=var(y, na.rm=T))Là encore, le mois de l’année ne semble guère influent vis à vis de ces caractéristiques, voir par exemple le tableau pour l’échéance à quatre jours.
d %>% filter(Echeance ==4) %>% mutate(y=log(V2/(T_std_ref_lis^2))) %>% group_by(Month) %>% summarise(moygeomvariance=exp(mean(y, na.rm=T)), VarLogVariance=var(y, na.rm=T))2.3.3 Relation entre variance de l’ensemble et écart de la température observée à la moyenne d’ensemble
On observe enfin un lien de l’écart entre observation et moyenne d’ensemble à la variabilité inter-membre. La figure ci-après illustre cette liaison (dans l’échelle des écart-types relatifs). Je reste pour le moment sans explication en ce qui concerne la valeur de la pente, proche de \(1\).
d %>% mutate(x=(V2/(T_std_ref_lis^2)),y=((Obs-Xbar)/T_std_ref_lis)^2) %>%
ggplot(aes(x=(x)^0.5,y=(y)^0.5))+geom_point(aes( color=as.factor(Month)))+geom_smooth(se=0, color='black')+geom_abline(intercept=0, slope=1, color='white')+facet_wrap(~Echeance)+labs(y="écart absolu prévision-temperature normalisé", x="std relatif", color="Mois")
On souhaite se servir de ce lien, même s’il s’avère assez flou, pour construire un modèle de prévision probabiliste reliant les statistiques résumées \(\bar{X}_t, V_t^2\) de l’ensemble et l’observation cible \(\theta_t\).
3 Le challenge
3.1 Objectifs
L’objectif de ce travail est:
d’utiliser la méthode cohérente du BFS pour construire une prévision probabiliste à partir du modèle échangeable normal à deux pivots,
de montrer ses performances sur les données de températures de Vouglans, à diverses échéances.
de comparer ses performances à celles des méthodes courantes EMOS et BMA, par exemple.
de montrer comment intégrer dans le même schéma les informations issues d’autres systèmes d’ensembles.
3.2 Les performances à battre
La figure ci-après annonce l’objectif à battre en terme de CRPS par rapport à la méthode d’ensemble du CEP.
label<-data.frame(
Echeance=c(8,5),
CRPS= c(2.05,1.1),
label=c("Climatologie","CEP")
)
d %>% filter(Year>2011) %>% group_by(Echeance) %>% mutate(crps_ens=crps(obs=Obs,pred = data.frame(Xbar,V2^0.5, na.rm=T))$crps,
crps_clim=crps(obs=Obs,pred = data.frame(T_moy_ref_lis,T_std_ref_lis, na.rm=T))$crps) %>% summarize(CRPS_ens=mean(crps_ens,na.rm=T),CRPS_clim=mean(crps_clim,na.rm=T)) %>%
ggplot(aes(x=Echeance))+geom_line(aes(y=CRPS_clim),color='red')+geom_line(aes(y=CRPS_ens))+geom_point(aes(y=CRPS_clim),color='red')+geom_point(aes(y=CRPS_ens))+labs(y="CRPS")+geom_label(data=label,aes(y=CRPS, label=label))
On cherchera également à vérifier la calibration, par exemple en testant le caractère uniforme de leur probability integral transforms.
d %>% filter(Year>2011) %>% mutate(pit_ens=crps(obs=Obs,pred = data.frame(Xbar,V2^0.5, na.rm=T))$pit,
pit_clim=crps(obs=Obs,pred = data.frame(T_moy_ref_lis,T_std_ref_lis, na.rm=T))$pit) %>%
ggplot()+geom_histogram(mapping=aes(x=pit_clim),stat="density", color='red')+
geom_histogram(mapping=aes(x=pit_ens),stat="density",position = "dodge", color="blue", alpha=0.3)+geom_hline(yintercept = 1)+facet_wrap(~Echeance)+labs(x="Probability Integral Transform")Ignoring unknown parameters: binwidth, bins, padIgnoring unknown parameters: binwidth, bins, pad
3.3 Diviser l’échantillon en apprentissage+validation
Pour tester la performance des diverses méthodes, on réalisera l’inférence sur un échantillon d’apprentissage (par exemple 2005-2008) et on validera sur l’échantillon restant (par exemple 2011-2015)
d %>% filter(Year < 2010, !is.na(Calendaire)) %>%
mutate(theta=(Obs-T_moy_ref_lis)/T_std_ref_lis,
xbar=(Xbar-T_moy_ref_lis)/T_std_ref_lis,
v2= V2/(T_std_ref_lis^2) ) %>%
filter(!is.na(theta),
!is.na(xbar),
!is.na(v2)) ->
learning_sample
d %>% filter(Year > 2010 ,!is.na(Calendaire)) %>%
mutate(theta=(Obs-T_moy_ref_lis)/T_std_ref_lis,
xbar=(Xbar-T_moy_ref_lis)/T_std_ref_lis,
v2= V2/(T_std_ref_lis^2) ) %>%
filter(!is.na(theta),
!is.na(xbar),
!is.na(v2)) ->
validation_sample4 Modèles d’ensemble
4.1 Modèle Gamma-Normal à deux pivots
Rappelons la suggestion de Jacques dans son document EnsembleSaison. Les membres s’appuient sur deux pivots latents \(Z_{1t}\) et \(Z_{2t}\), de telle sorte que le modèle échangeable s’écrit:
\[{X_{ts}} = \alpha + \beta {Z_{1t}} + \lambda \sigma {Z_{2t}}^{ - O.5}{\varepsilon _{ts}}\]
où \(\varepsilon _{ts}\sim N(0,1)\) tandis que \(\sigma\) peut être pris égal à \(1\). Les pivots latents sont munis de priors \[\begin{gathered} {Z_{2,t}} \sim Gamma(g,1) \hfill \\ {Z_{1t}}|{Z_{2t}} \sim N(0,{\sigma ^2}{Z_{2t}}^{ - 1}) \hfill \\ \end{gathered}\]
Travaillant avec les statistiques exhausives, on écrira \[\begin{gathered} {{\bar X}_t} \sim N(\alpha + \beta {Z_{1t}},\frac{{{\lambda ^2}{\sigma ^2}}}{{S{Z_{2t}}}}) \hfill \\ \sum\limits_{s = 1}^S {{{({X_{ts}} - {{\bar X}_t})}^2}} \sim \frac{{{\lambda ^2}{\sigma ^2}}}{{{Z_{2t}}}}{\chi ^2}(S - 1) \hfill \\ \end{gathered}\]
Soit encore, en posant \(V_{t}^{2}=\frac{{\sum\limits_{s = 1}^S {{{({X_{ts}} - {{\bar X}_t})}^2}} }}{{S - 1}}\), le modèle d’échantillonnage à \(t\), par orthogonalité des deux statistiques exhaustives moyenne et variances empiriques, se résume par les deux équations: \[\begin{gathered} {{\bar X}_t} \sim N(\alpha + \beta {Z_{1t}},\frac{{{\lambda ^2}{\sigma ^2}}}{{S{Z_{2t}}}}) \hfill \\ {V_{t}^{2}} \sim Gamma\left( {\frac{{S - 1}}{2},\frac{{{Z_{2t}}(S - 1)}}{{2{\lambda ^2}{\sigma ^2}}}} \right) \hfill \\ \end{gathered}\]
L’inférence des coefficients \(\alpha ,\beta ,\lambda , g\) sera menée sur l’ échantillon d’apprentissage. Prenons par exemple l’écheance à 4 jours des années 2005-2008.
learning_sample %>% filter( Echeance == 4) %>%
dplyr::select(theta,xbar,v2) -> apprentissageOn procède à l’estimation bayésienne de \(\alpha ,\beta ,\lambda , g\) grâce au logiciel d’inférence bayésienne Jags. (A noter qu’on pourrait également utiliser STAN.)
model_string <- "
model{
a<-(S-1)/2
for (t in 1:N){
m[t] <- alpha+beta*Z1[t]
preci[t]<- Z2[t]/(lambda2)
xbar[t] ~ dnorm(m[t], preci[t])
b[t]<-preci[t]*a
v2[t] ~ dgamma(a,b[t])
# latentes
Z2[t] ~ dgamma(g,1)
Z1[t] ~ dnorm(0, 1)
}
alpha ~ dunif(-10,10)
beta ~ dunif(0,10)
lambda2 ~ dunif(0,10)
g ~ dunif(0.5,10)
}
"
params=c("alpha","beta","lambda2","g")
data <- list(xbar=apprentissage$xbar,
v2=apprentissage$v2,
N=length(apprentissage$xbar),
S=50)
temp<-jags(data = data, model.file = textConnection(model_string),
parameters.to.save = params,n.chains = 3,
n.burnin = 5000,n.iter = 10000)Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 2908
Unobserved stochastic nodes: 2912
Total graph size: 11646
Initializing model
|
| | 0%
|
|++ | 4%
|
|++++ | 8%
|
|++++++ | 12%
|
|++++++++ | 16%
|
|++++++++++ | 20%
|
|++++++++++++ | 24%
|
|++++++++++++++ | 28%
|
|++++++++++++++++ | 32%
|
|++++++++++++++++++ | 36%
|
|++++++++++++++++++++ | 40%
|
|++++++++++++++++++++++ | 44%
|
|++++++++++++++++++++++++ | 48%
|
|++++++++++++++++++++++++++ | 52%
|
|++++++++++++++++++++++++++++ | 56%
|
|++++++++++++++++++++++++++++++ | 60%
|
|++++++++++++++++++++++++++++++++ | 64%
|
|++++++++++++++++++++++++++++++++++ | 68%
|
|++++++++++++++++++++++++++++++++++++ | 72%
|
|++++++++++++++++++++++++++++++++++++++ | 76%
|
|++++++++++++++++++++++++++++++++++++++++ | 80%
|
|++++++++++++++++++++++++++++++++++++++++++ | 84%
|
|++++++++++++++++++++++++++++++++++++++++++++ | 88%
|
|++++++++++++++++++++++++++++++++++++++++++++++ | 92%
|
|++++++++++++++++++++++++++++++++++++++++++++++++ | 96%
|
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
|
| | 0%
|
|** | 4%
|
|**** | 8%
|
|****** | 12%
|
|******** | 16%
|
|********** | 20%
|
|************ | 24%
|
|************** | 28%
|
|**************** | 32%
|
|****************** | 36%
|
|******************** | 40%
|
|********************** | 44%
|
|************************ | 48%
|
|************************** | 52%
|
|**************************** | 56%
|
|****************************** | 60%
|
|******************************** | 64%
|
|********************************** | 68%
|
|************************************ | 72%
|
|************************************** | 76%
|
|**************************************** | 80%
|
|****************************************** | 84%
|
|******************************************** | 88%
|
|********************************************** | 92%
|
|************************************************ | 96%
|
|**************************************************| 100%temp$BUGSoutput$sims.matrix %>%
as.data.frame() %>%
dplyr::select(alpha,beta,lambda2,g) ->alphabetalambda2g
alphabetalambda2g %>%
gather(param, value) %>%
ggplot(aes_string(x="value")) +
geom_density(alpha=0.5) +
facet_grid(param~.,scales = "free")
L’intercept \(\alpha\) est possiblement nul, la pente \(\beta\) restant inférieure à \(1\). La connaissance a posteriori de \(g\) est assez incertaine; \(\lambda^2\) et \(g\) sont fortement corrélés.
ggpairs(data = alphabetalambda2g)
knitr::kable(rbind(mean=apply(X = alphabetalambda2g,2,mean),std=apply(X = alphabetalambda2g,2,var)^0.5,apply(X = alphabetalambda2g,2, quantile, probs=c(0.05,0.25,0.5, 0.75, 0.95))), dig=3)Registered S3 methods overwritten by 'htmltools':
method from
print.html tools:rstudio
print.shiny.tag tools:rstudio
print.shiny.tag.list tools:rstudio
| alpha | beta | lambda2 | g | |
|---|---|---|---|---|
| mean | 0.023 | 0.904 | 0.097 | 1.778 |
| std | 0.025 | 0.018 | 0.004 | 0.063 |
| 5% | -0.018 | 0.875 | 0.090 | 1.676 |
| 25% | 0.006 | 0.891 | 0.094 | 1.735 |
| 50% | 0.023 | 0.903 | 0.097 | 1.777 |
| 75% | 0.039 | 0.916 | 0.100 | 1.821 |
| 95% | 0.064 | 0.933 | 0.104 | 1.880 |
knitr::kable(cor(alphabetalambda2g),dig=3)| alpha | beta | lambda2 | g | |
|---|---|---|---|---|
| alpha | 1.000 | -0.108 | -0.022 | -0.016 |
| beta | -0.108 | 1.000 | 0.087 | 0.082 |
| lambda2 | -0.022 | 0.087 | 1.000 | 0.871 |
| g | -0.016 | 0.082 | 0.871 | 1.000 |
4.2 Analyse marginale du modèle à deux pivots
La vraisemblance à \(t\) s’écrit: \[\begin{gathered} {L_t} = [{{\bar X}_t},{V_t}^2|{Z_{1t}},{Z_2}] \hfill \\ {L_t} \propto {Z_2}^{\frac{1}{2}}\exp - \frac{{{Z_{2t}}S}}{{2{\lambda ^2}}}\left\{ {{{\left( {{{\bar X}_t} - \alpha - \beta {Z_{1t}}} \right)}^2}} \right\} \times {Z_{2t}}^{\frac{S}{2} - \frac{1}{2}}\exp - \frac{{{Z_{2t}}}}{{2{\lambda ^2}}}\left\{ {\left( {S - 1} \right){V_t}^2} \right\} \hfill \\ \end{gathered}\]
On va travailler comme Krzysztofowicz en négligeant les incertitudes des hyper-paramètres du modèle marginal.
alpha<-mean(alphabetalambda2g[,"alpha"])
beta<-mean(alphabetalambda2g[,"beta"])
lambda2<-mean(alphabetalambda2g[,"lambda2"])
g<-mean(alphabetalambda2g[,"g"])
S<-50
Repet <-1000
data.frame(Z2=rgamma(Repet,g,1), Z1=rnorm(Repet)) %>%
mutate(xbar=alpha+beta*Z1+sqrt(lambda2/Z2/S)*rnorm(Repet),
v2=lambda2/Z2/(S-1)*rchisq(Repet,df=S-1)) %>% filter(v2<2) %>%
dplyr::select(xbar,v2) %>% gather(key,value) %>% ggplot(aes(x=value,color=key))+
geom_histogram(data= apprentissage %>% dplyr::select(xbar,v2) %>% filter(v2<2) %>% gather(key,value), stat = "density")+
geom_freqpoly(stat = "density", color='black')+facet_wrap(~key,scales="free")Ignoring unknown parameters: binwidth, bins, pad
Sur cette comparaison estimation empirique et modèle, on voit que les lois marginales se calent raisonnablement avec notre modèle gamma-normal echangeable, malgré que nou sn’ayons pas tenu compte des incertitudes à propos des hyper-paramètres (\(\alpha,\beta, \gamma^2, g\)) du modèle marginal.
4.3 Analyse préquentielle du modèle à deux pivots
4.3.1 Expression de la conditionnelle des deux variables pivots latentes
La vraisemblance complète quant à elle est légèrement plus compliquée: \[\begin{gathered} C{L_t} = [{{\bar X}_t},{V_t}^2,{Z_{1t}},{Z_2}] = [{{\bar X}_t},{V_t}^2|{Z_{1t}},{Z_{2t}}] \times [{Z_{1t}},{Z_{2t}}] \hfill \\ C{L_t} \propto {L_t} \times \left( {{Z_2}^{\frac{1}{2}}\exp - \frac{{{Z_{2t}}}}{2}Z_{1t}^2} \right) \times \left( {Z_{2t}^{g - 1}\exp - {Z_{2t}}} \right) \hfill \\ \end{gathered}\]
En développant: \[\begin{gathered} C{L_t} \propto \left( {{Z_2}^{\frac{1}{2}}\exp - \frac{{{Z_{2t}}}}{2}\left\{ {Z_{1t}^2 + \frac{{S{\beta ^2}}}{{{\lambda ^2}}}{{\left( {{Z_{1t}} - \frac{{{{\bar X}_t} - \alpha }}{\beta }} \right)}^2}} \right\}} \right) \hfill \\ \times \left( {Z_{2t}^{g + \frac{S}{2} - 1}\exp - {Z_{2t}}(1 + \frac{{\left( {S - 1} \right){V_t}^2}}{{2{\lambda ^2}}})} \right) \hfill \\ \end{gathered}\]
puis en regroupant:
\[\begin{gathered} Z_{1t}^2 + \frac{{S{\beta ^2}}}{{{\lambda ^2}}}{\left( {{Z_{1t}} - \frac{{{{\bar X}_t} - \alpha }}{\beta }} \right)^2} = \left( {\frac{{S{\beta ^2} + {\lambda ^2}}}{{{\lambda ^2}}}} \right){\left( {Z_{1t}^{} - \left( {\frac{{{{\bar X}_t} - \alpha }}{\beta }} \right)\frac{{S{\beta ^2}}}{{S{\beta ^2} + {\lambda ^2}}}} \right)^2} \hfill \\ + \left\{ {\left( {\frac{{S{{\left( {{{\bar X}_t} - \alpha } \right)}^2}}}{{S{\beta ^2} + {\lambda ^2}}}} \right)} \right\} \hfill \\ \end{gathered}\]
la conjugaison adhoc permet d’obtenir les conditionnelles a posteriori
\[\begin{gathered} {Z_{2t}}|{{\bar X}_t},{V_t}^2 \sim Gamma(g + \frac{S}{2},1 + \frac{{\left( {S - 1} \right){V_t}^2}}{{2{\lambda ^2}}} + \left( {\frac{{S{{\left( {{{\bar X}_t} - \alpha } \right)}^2}}}{{2\left( {S{\beta ^2} + {\lambda ^2}} \right)}}} \right)) \hfill \\ {Z_{1t}}|{Z_{2t}},{{\bar X}_t},{V_t}^2 \sim N(\left( {\frac{{{{\bar X}_t} - \alpha }}{\beta }} \right)\frac{{S{\beta ^2}}}{{S{\beta ^2} + {\lambda ^2}}},{\left\{ {\left( {\frac{{S{\beta ^2} + {\lambda ^2}}}{{{\lambda ^2}}}} \right){Z_{2t}}} \right\}^{ - 1}}) \hfill \\ \end{gathered}\]
4.3.1.1 Nota:
Je pense qu’ il y a une erreur typographique dans le papier de Jacques en page 8 pour le calcul de l’espérance de \(Z_1\).
4.3.2 Apprentissage d’un lien des deux variables pivots latentes avec la variable cible
Pour chaque situation de la base d’apprentissage, nous générons un \(Z_1(t)\) que nous mettons en regard avec la cible.
La suggestion de Jacques est de repasser en marginale et de travailler avec les QNT:
la cible qui a été centrée réduite est déjà N(0,1).
\(Z_1\) est marginalement Student à \(2g\) degrés de liberté. Quant à \(g^{0.5} \times Z_1\), il est Student unitaire à \(2g\) degrés de liberté.
lm(QNT_Z1~theta-1, data=apprentissage)$coef
theta
0.8803137 Sur l’échantillon d’apprentissage, on trouve un coefficient de corrélation \(\hat{\rho}=0.8803\) entre le QNT de la cible et le QNT d’une valeur tirée au hasard dans la loi conditionnelle du premier pivot sachant les deux statistiques résumant l’ensemble échangeable.
4.4 Analyse prédictive grâce au modèle à deux pivots
Considérons maintenant l’échantillon de validation ou l’on va appliquer le calcul prédictif:
\([\theta|X]=\int_{QNT} [\theta|QNT] \times [QNT|X] \times dQNT\)
Et si l’on dispose d’un prior informatif \([\theta]\) pour la cible:
\[[\theta|X]=\int_{QNT} \frac{[QNT|\theta]\times [\theta]}{[QNT]} \times [QNT|X] \times dQNT\] Si on utilise la conjugaison normale avec un prior \(N(m_t,s_t^2)\) pour \([\theta_t]\), \([\theta|QNT]\) est \(N(m'_t,s'_t^2)\) avec
\(\frac{1}{s'^2}=\frac{1}{s^2}+\frac{\rho^2}{1-\rho^2}\)
\(\frac{m'}{s'^2}=\frac{m}{s^2}+\frac{QNT}{\rho}\frac{\rho^2}{1-\rho^2}\)
pour \(m=0, s=1\), c’est à dire pour un prior climatologique pour produire une prévision probabiliste calibrée, on retrouve les résultats de la régression \([\theta|QNT]\) est \(N(\rho \times QNT,1-\rho^2)\)
validation_sample %>% filter(Echeance==4) %>% mutate(a_Z2=g+S/2, b_Z2=1+(S-1)*v2/2/lambda2+S*((xbar-alpha)^2)/2/(S*beta^2+lambda2),
m_Z1= (xbar-alpha)*S*beta /(S*beta^2+lambda2) ,
precision_Z1=(S*beta^2+lambda2)/lambda2) ->validation 4.5 EMOS
Je veux comparer à un modèle bayésien qui conditionne les valeurs de l’ensemble sur la température vraie que j’ai appellé \(\theta\) car c’ est la notation conventionnelle de l’état de la nature d’un problème statistique. Je propose : \[ \theta \sim N(0,1) \\ \bar{X}\vert (V^2,\theta) \sim N(a\times \theta +b , c+ d\times V^2) \\ log(V^2) \sim N(g,s^2) \]
L’inférence des coefficients \(a,b,c,d,g,s^2\)) sera menée sur un échantillon d’apprentissage.
La loi prédictive (c-à-d conditionnelle sachant les résumés de l’ensemble \(\bar{X}\) et \(V^2\)) est ici “immédiate”, par conjugaison normale:
\[\theta \vert (\bar{X},V^2,a,b,c,d) \sim N\left\{ {\left( {\frac{{\bar X}}{a} - b} \right)\frac{{{a^2}}}{{{a^2} + c + d{V^2}}},\frac{{c + d{V^2}}}{{{a^2} + c + d{V^2}}}} \right\}\]
Ses performances seront analysées sur un échantillon de validation. Prenons par exemple l’écheance à 4 jours des années 2005-2008.
d %>% filter(Year < 2010, Echeance == 4, !is.na(Calendaire)) %>%
mutate(theta=(Obs-T_moy_ref_lis)/T_std_ref_lis,
xbar=(Xbar-T_moy_ref_lis)/T_std_ref_lis,
v2= V2/(T_std_ref_lis^2) ) %>%
dplyr::select(theta,xbar,v2) %>% filter(!is.na(theta),
!is.na(xbar),
!is.na(v2)) ->
learning_sampleOn va réaliser l’inférence bayésienne des coefficients \(a,b,c,d\)
model_string <- "
model{
for (i in 1:N){
m[i] <- a*theta[i]+b
preci[i]<- 1/(c+d*v2[i])
xbar[i] ~ dnorm(m[i], preci[i])
}
a ~ dunif(-10,10)
b ~ dunif(-10,10)
c ~ dunif(0,10)
d ~ dunif(0,10)
}
"
params=c("a","b","c","d")
data <- list(theta=learning_sample$theta,
xbar=learning_sample$xbar,
v2=learning_sample$v2,
N=length(learning_sample$theta))
temp<-jags(data = data, model.file = textConnection(model_string),
parameters.to.save = params,n.chains = 3,
n.burnin = 500,n.iter = 1000)temp$BUGSoutput$sims.matrix %>%
as.data.frame() %>% dplyr::select(a,b,c,d) ->abcd
abcd %>% gather(param, value) %>%
ggplot(aes_string(x="value")) +
geom_density(alpha=0.5) +
facet_grid(param~.,scales = "free_x")
Les coefficients \(a,b,c\) sont assez précisément estimés sauf la valeur de \(d\) qui est plus incertaine, tandis que les coefficients \(c\) et \(d\) sont anti-corrélés.
knitr::kable(rbind(mean=apply(X = abcd,2,mean),std=apply(X = abcd,2,var)^0.5,apply(X = abcd,2, quantile, probs=c(0.05,0.25,0.5, 0.75, 0.95))), dig=3)
knitr::kable(cor(abcd),dig=3)4.6 Predictive bayésienne
Sur un échantillon de validation (ici 2011-2015), on va étudier le comportement du prédicteur bayésien. \[\theta \vert (\bar{X},V^2,a,b,c,d) \sim N\left\{ {\left( {\frac{{\bar X}-b}{a}} \right)\frac{{{a^2}}}{{{a^2} + c + d{V^2}}},\frac{{c + d{V^2}}}{{{a^2} + c + d{V^2}}}} \right\}\]
d %>% filter(Year > 2010, Echeance == 4, !is.na(Calendaire)) %>%
mutate(theta=(Obs-T_moy_ref_lis)/T_std_ref_lis,
xbar=(Xbar-T_moy_ref_lis)/T_std_ref_lis,
v2= V2/(T_std_ref_lis^2) ) %>%
filter(!is.na(theta),!is.na(xbar),!is.na(v2)) ->
validation_sample
indices<-sample(x = 1:dim(abcd)[1], size=50, rep=F)
A=abcd[indices,1]
B=abcd[indices,2]
m=(outer(validation_sample$xbar,A,'/')-B/A)
C=abcd[indices,3]
D=abcd[indices,4]
s=outer(validation_sample$v2,D,'*')+C
m=m*A*A/(A^2+s)
s=(s/s+A*A)^0.5
YPRED=matrix(
rnorm(n = prod(dim(s)), mean = m, sd = s),
nr=dim(s)[1],nc=dim(s)[2] )*
validation_sample$T_std_ref_lis+validation_sample$T_moy_ref_lis
moyPRED=apply(YPRED, 1, mean)
stdPRED=apply(YPRED, 1, sd)
rm(YPRED)